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Abstract 

We performed numerical simulations of the obliquity evolution of Mars during the Noachian era, at which time the 
giant planets were on drastically different orbits than today. For the preferred primordial configuration of the planets 
we find that there are two large zones where the Martian obliquity is stable and oscillates with an amplitude lower than 
20°. These zones occur at obliquities below 30° and above 60°; intermediate values show either resonant or chaotic 
behaviour depending on the primordial orbits of the terrestrial planets. 
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1. Introduction 

The current obliquity of Mars is chaotic with a Lyapunov time scale of 5 Myr (Laskar & Robutel, 1993; Touma & 
Wisdom, 1993). This has drastic consequences for the Martian climate. At present, the obliquity is a modest 25°, the 
average polar insolation is lower than the equatorial one (Ward, 1974) and ice caps are observed at the poles of the 
planet. However, during phases where the obliquity is high, e > 40°, large quantities of polar ice sublimate and could 
be transported to the tropical regions (Jakosky & Carr, 1985; Jakosky et al., 1995). On shorter time scales the Mar- 
tian obliquity suffers large-amplitude oscillations of approximately 22° around a mean value of 25° with a period of 
125 kyr (Ward, 1974), caused by the perturbations on its orbit from the other terrestrial planets. These large-amplitude 
oscillations drive large-scale variations in its seasonal cycles of water, carbon dioxide and dust (e.g. Kieffer & Zent, 
1992). In addition, the large obliquity changes cause frequent ice ages on Mars, during which there is a build-up of 
polar ice towards the equator (Head et al, 2003). The planet is currently in an interglacial period having suffered an 
ice age approximately 1 Myr ago when the obliquity reached over 30° (Head et al, 2003). Over the past 10 Myr the 
mean obliquity has decreased from approximately 35° to its current value of 25° (Touma & Wisdom, 1993; Laskar et 
al. 2002), indicating that past ice ages may have lasted much longer. 

However, the conditions on Mars were very different shortly after its formation compared to what they are to- 
day. During the Noachian period, which lasted from Mars' formation until 3.5 Gyr ago, there is evidence for a short 
existence of liquid water during its later stages (e.g. Carr, 1996; Christensen et al., 2001), a temporary magnetic 
field (Lillis et al., 2008) and possibly a thick atmosphere with the temperature at the surface exceeding 273 K. The 
existence of ubiquitous valley networks has led to the idea that fluvial erosion was an important process on Noachian 
Mars (Carr, 1996). However, Christensen et al. (2001) argue that the duration of liquid water flow must have been 
brief because of the existence of olivine and lack of carbonates at the surface.By the end of the Noachian, the bulk of 
the valley networks had formed and the erosion rate was steeply decUning (Carr, 1996). In addition, the magnetic field 
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might have protected the atmosphere from pick-up-ion sputtering and hydrodynamic colhsions, and therefore knowl- 
edge of its evolution is crucial for understanding the Noachian climate history. Isotopic evidence in Martian meteorite 
ALH84001 suggests that the atmosphere was largely unfractionated near the end of the Noachian (Marti et al., 2001; 
Mathew & Marti, 2001). The geologic and isotopic information taken together suggests a relatively rapid loss of 
atmosphere as Mars entered the Hesperian era - which lasted from 3.5 Gyr ago to 1.8 Gyr ago. Loss of a global mag- 
netic field has been considered as an attractive mechanism, but this probably happened much earlier (Lillis et al, 2008). 

The end of the Noachian age corresponds approximately with the beginning of the Lunar Cataclysm or Late 
Heavy Bombardment (LHB). This episode of intense cratering on the Moon (Tera et al., 1974), and the other terres- 
trial planets, occurring some 3.8 Gyr ago, formed some of the basins on the Moon, notably Imbrium and Orientale. 
The cratering intensity was dramatically higher than what is predicted from the remnants of planet formation (Bottke 
et al., 2007), suggesting a major re-ordering event in the solar system. Indeed, the LHB is beUeved to be the last 
major event that sculpted the solar system (Strom et al., 2005). The LHB is also thought to be linked with a dynam- 
ical instability in the outer solar system, which resulted in a reshuffling of giant planet orbits (Gomes et al., 2005). 
This dynamical instability also affected the orbits of the terrestrial planets and probably created their current orbital 
structure (Brasser et al., 2009). These models favour a fast migration of Jupiter and Saturn, resulting in a short-Uved, 
high-eccentricity phase after a long spell of circular orbits (Brasser et al., 2009). Recent results have uncovered that 
the orbits of the giant planets were likely less dynamically excited (lower eccentricities and inclinations), and that 
they were much closer together (Morbidelli et al., 2007). However, for this study we shall not limit ourselves to this 
one particular case, but explore a wider range of possible configurations. Unfortunately, the orbits of the terrestrial 
planets before the LHB are unknown, but there is a range of possible solutions. Therefore we create a set of different 
terrestrial planet systems with a wide range of initial conditions, and determine the influence of their perturbations 
on the obliquity evolution of Mars using numerical simulations. The simulations are inspected for cases where the 
obliquity could be stable on long time scales. 

This paper proceeds as follows. Section 2 contains a short presentation of the equations of motion of the obliquity 
of Mars, and the typical values of the quantities that enter these equations. Section 3 contains a description of our 
numerical methods. In Section 4 we present our results, and the last section contains our conclusions. 



2. Evolution of the obliquity with planetary perturbations 

The evolution of the obliquity of a planet subject to perturbations from other planets has been studied by various 
authors (e.g. Colombo, 1966; Ward, 1974; Kinoshita, 1977; Laskaret al., 1993; Ward & Hamilton, 2004). Here we 
give a short summary based on the elegant method of Neron de Surgy & Laskar (1997). 

Mars' figure exhibits precession of the equinoxes, similar to Earth. The precession is caused by Mars' figure being 
slightly flattened, so that there is a normal component to the force exerted on Mars by the Sun. This component of 
the force induces a torque on Mars that drives the precession and which acts in the direction opposite to Mars' orbital 
motion. The Hamiltonian describing the torque on Mars is (Neron de Surgy & Laskar, 1997) 

(1) 

2C 2L ^ ^ 

where L = Cw, is the rotational angular momentum of Mars, C is Mars' largest moment of inertia, co,- is its rotation 
rate, X - L cos s with e being the obliquity. The precession of the equinox, i/r, is governed by the parameter a as 
- aX where 

2a)r C 



Here J2 is Mars' quadrupole moment, n - ^IGMqIo^ is Mars' orbital mean motion, Rm ^^ Mars' radius, M is its mass, 
a is the semi-major axis and e is the eccentricity. For simplicity we will set L = 1 and thus X - cos e. The most recent 
data for the figure, obliquity, equinox and orbit of Mars are taken from Folkner et al. (1997), with e = 25.189517°, 



2 



(A = 35.43777°, = 350.891985217day, C/M/?^ = 0.3662, Jx = 1.96045 x 10-^ a = 1.523698 AU and 
e - 0.0933428, from which we compute the current value a = 8.373 "/yi and thus acos s - 7.576 "/yr. 

However, the orbit of Mars is not static but undergoes periodic changes induced by perturbations from the other 
planets. When taking these perturbations into account, the ecliptic is no longer an inertial plane and thus the kinetic 
energy of its forcing, £, has to be added to the Hamiltonian (Kinoshita, 1977). This is given by 



6 = 2r{t)X - Vl - X^A(t) sin lA + B(t) cos i/r], (3) 

where 

r(t)^pq-qp; A(t) ^ 2— ; B{t) ^ 2—— (4) 

with ( - q + \p - sin(//2) exp(iQ) and i^ = -1. Here ; is the inclination of the orbit of Mars with respect to the 
invariable plane of the solar system and is its longitude of the ascending node. The new Hamiltonian is then 
TC = 'H + £, and Hamilton's equations are 



X = = - Vl -X2[A(f) COS lA-BCO sin i/r]; = = aZ-Z(l -XV/2[A(f) sintA + B(f) cos i/r] -2r(f). (5) 

oiA aX 

Here A(f), B{t) and r(f) are the forcing terms from the planetary perturbations on the orbit of Mars. Equations (|5) 
are singular when X - \ i.e. when e = 0, so for numerical integrations we follow Laskar et al. (1993) and use 
= ^ + i77 = sin e exp(ii/'), which yields 

^ = A{t) ^1 - ^2 _ ^2 _ ^1 _ ^2 _ ^2 _ 2r(f)]; 77 = -fi(0 ^l-e-ri^+ ^^-e-V^- 2r(f)]. (6) 

Equations (l6]i contain no singularities. Now that we have outlined the basic equations of motion of the spin axis, we 
turn to our numerical methods. 



3. Numerical Methods 

Unlike Laskar et al. (2004), our aim is not to obtain an accurate prediction of the future evolution of the obliquity 
of Mars, but rather to test its stability as a function of varying initial conditions: different orbits of the terrestrial 
planets and a range of original obliquities. Thus we integrate the equations of motion of the obliquity of Mars (|6]l 
from existing simulations of the solar system. This is described below. 

• We perform simulations of the solar system with various initial conditions of the giant planets and the terrestrial 
planets (see next section) for 265 Myr with a time step of 0.01 yr (3.6525 days). The data is output every 
1000 yr We use the 2nd order MVS integrator of Laskar & Robutel (2001), implemented (Broz et al., 2005) 
in the SWIFT integration package (Levison & Duncan, 1994). The code includes digital filtering of the orbital 
elements using Kaiser windows (Kaiser, 1966), based on the method of Quinn et al. (1991). Any signal with a 
period shorter than 667 yr is suppressed. To mimic the effect of general relativity we added a correction term as 
described in Nobili & Will (1986). 

• From the simulation the quantities a, A{f), B{f) and T{f) are calculated using the filtered elements. Numerical 
forward differentiation with four data points was used to compute p and q. 

• The quantities a, A(f), B{f) and Y{t) are used to integrate eqs. (|6) with the Bulirsch-Stoer method (Bulirsch 
& Stoer, 1966) with a time step of 1000 years. For the intermediate steps that the integrator performs before 
extrapolation to zero step size the corresponding intermediate values of a, A(f), B{t) and r(f) are computed using 
linear interpolation. The input are the original values of e and i/r. The output is Fourier analysed according to the 
method of Laskar (1988). The maximum and minimum values of e and the averaged proper rate of precession, 
^, are recorded. 
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Figure I : Top left panel shows maximum and minimum values of obliquity evolution (vertical axis) taken over 256 Myr as a function of the original 
value (horizontal axis). Top right panel shows the averaged proper precession rate, i// (vertical axis) vs original obliquity. The bottom panel shows 
a spectrum of A(t) + \B(t), taken over 64 Myr The vertical axis con'esponds to the log of the amplitude of each term while the horizontal axis 
corresponds to the frequency in arcsec per year. This figure is similar to those presented in Laskar & Robutel (1993). 



The last two steps are repeated for different values of original e. Generally e was varied between and 90° with 
steps of 0.2°, while the original value of was kept at its current value. All steps were repeated for the different 
configurations of the terrestrial planets. 



4. Results 

4.1. Current solar system 

In order to be able to distinguish the obliquity evolution of Mars during the Noachian from the current one, we 
give a short summary of the latter For more details see Laskar & Robutel (1993), Touma & Wisdom (1993) and 
Laskar et al., (2004). 

Fig. [U depicts current obliquity evolution, and should be compared directly with those presented in Laskar & 
Robutel (1993). The top left panel of Fig. [T] depicts the minimum and maximum values that the obliquity reaches 
during 256 Myr (vertical axis) as a function of initial obliquity, bq (horizontal axis). The top right panel depicts the 
average proper precession rate of the equinoxes, i/^, on the vertical axis as a function of sq. There are three regions of 
motion. The first is for initial obliquities < 10°, where it oscillates between and 20° and the precession decreases 
with the obliquity as i/' = a cose. The second region is between 10° and 60°. Here the minimum and maximum 
obliquities are 10° and 60° respectively and ranges from 5 "/yr to 8 "/yr without a trend with e. This is the chaotic 
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region identified by Laskar & Robutel (1993) and Touma & Wisdom (1993). On Gyr time scales, the mean obliquity 
of Mars will vary between and 60°. Above 60° the motion becomes regular. The bottom panel depicts the Fourier 
spectrum of A{t) + iB{t), similar to Laskar & Robutel (1993), but taken over 64 Myr The vertical axis lists the log 
of the amplitude of each term while the horizontal axis corresponds to the frequency in arcsec per year. The forcing 
peaks from the solar system's nodal eigenfrequencies si . . . s^, corresponding to the proper regression frequencies of 
the nodes of Mercury (si) to Neptune (ig), (Brouwer & van Woerkom, 1950), are indicated. 

It is important to distinguish the two sources of chaos in the current evolution of Mars' obliquity. The first is 
caused by a resonance overlap between the precession of Mars' figure and two eigenfrequencies of the solar system. 
The other is a secular resonance between the orbits of Mercury, Venus and Jupiter, which indirectly affects the orbit 
of Mars. We shall examine each of these below. 

Usually the relation ij/ - a cose is a good approximation for the precession rate with increasing obliquity. For 
Mars, however, as the obliquity increases there are two values where ij/ = i.e. a secular spin-orbit resonance occurs. 
The first resonance is with ^2 at S2 - arccos(i2/Q;) = 32.3°. The second is with si at ei - arccos(ii/Q') = 48.0°. 
The width of these two resonances is directly proportional to the forcing on Mars' orbit from Mercury and Venus. 
In the current solar system the forcing from both planets is so large that these resonances overlap, which occurs at 
s ~ 40°. Resonance overlap is a source of chaos (Chirikov, 1979), which accounts for the chaotic region between 
initial obliquities of 30° and 60°. However, there is a second source of chaos that affects the obliquity, which arises 
from chaos in the orbit of Mars. 

Laskar (1990) identified two sources of chaos in the orbits of the terrestrial planets. The one that affects Mars' 
obliquity has the argument cr = {gi - gs) - {s[ - S2), which is currently in libration with period 10 Myr (Laskar, 
1990). Here gi . . .gs. are the eigenfrequencies of the longitudes of pericentre, corresponding to the proper precession 
frequencies of the perihelia of Mercury (gi) to Neptune (g^) (Brouwer & Van Woerkom, 1950). The argument cr 
stays in libration for the next 200 Myr, but its amplitude changes, implying the motion is near a separatrix. Laskar 
(1990) showed that the components associated with si and S2 in Mars' orbit are accompanied by a number of smaller 
components of similar frequency. A multiplet of components can also be viewed as a single component with varying 
frequency and amplitude, with the frequency of the multiplet being comparable to the frequency spread of the multi- 
plet. These multiplets are caused by the presence of the secular resonance {gi - gs) - (ii - ^2). Thus in the Fourier 
spectrum of Mars' orbit the lines associated with si and S2 are not clean delta functions but exhibit a near-Gaussian 
profile i.e. they have significant sidebands. This is clearly visible in fig.[T]and the profiles overlap. These mutiplets 
effect the motion of the mean obliquity because ip can resonate with all these additional frequencies in the multiplets. 
The mean obliquity pops in and out of resonance because the location of the separatrix varies with time. The crossing 
of the separatrix occurs at essentially random phase which causes the obliquity to increase or decrease at random, 
and thus the mean obliquity exhibits something akin to a random walk. The secular resonance acts when the mean 
obliquity < 30°. This implies that the mean obliquity is chaotic from 0° to 60° on Gyr time scales (Laskar & Robutel, 
1993; Laskar et al., 2004). 

From the above it appears there are two ways the obliquity of Mars can be stabilised. The first is to reduce the 
forcing from Mercury and Venus on Mars so that the two resonances no longer overlap. This can be achieved by 
reducing the eccentricities and inclinations of the terrestrial planets. There is indication that the terrestrial planets 
were situated on dynamically colder orbits before the LHB (Brasser et al, 2009). The second method is to change the 
precession frequencies of either Mercury, Venus or Jupiter and break the resonance (gi - ^5) - (ii - 52)- There are no 
indications that the terrestrial planets had their orbits changed significantly, but there is plenty of evidence that Jupiter 
had migrated (e.g. Fernandez & Ip, 1984; Malhotra, 1993; Tsiganis et al., 2005; Morbidelh et al., 2007, 2009, 2010; 
Brasser et al., 2009). Displacing Jupiter would change the value of ^5 and if increased enough it will break the secular 
resonance. The effect of removing this resonance on the obliquity of Mars, while keeping everything else the same, is 
discussed in the next subsection. 
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Figure 2: Similar to Fig.^ but now the giant planets are on the configuration of Malhotra (1993, 1995). The top-left panel depicts the maximum 
and minimum obliquity as a function of original obliquity. The top-light panel shows the averaged proper precession frequency. Both two panels 
are taken over 256 Myr The bottom panel shows a spectrum of A(t) + iB(t), taken over 64 Myr Notice that the chaotic motion is confined to values 
of the initial obliquity between 35° and 60°. 



4.2. Compact configuration of the giant planets 

Here we investigate the effect of a different configuration of the giant planets on the obliquity evolution of Mars. 
It was almost certain that the giant planets had different orbits before the LHB, so it is essential to study the stability 
of Mars' obliquity using other configurations of the giant planets. We look for a system that increases ^5 so that the 
resonance (g] - gs) - (si - S2) is not in libration. 

The value of ^5 depends mostly on the separation between Jupiter and Saturn i.e. on their period ratio Ps /Pj- 
Brasser et al. (2009) showed that when Ps /Pj > 2, decreases sharply as Ps /Pj increases: from approximately 
20"/yr when PsIPj - 2.03 to its current value of 4.26"/yr when Ps /Pj = 2.48. As Ps /Pj increases and ^5 decreases, 
it crosses the values of §2 = 7.45"/yr and ^1 = 5.56'7yr when Ps /Py ~ 2.1 andPs/Pj ~ 2.3 respectively. In order 
not to destabilise the terrestrial system we need gs > g2 and thus Ps /Pj < 2.1. There are several well-studied models 
in the literature that prefer an initial period ratio between Jupiter and Saturn PsIPj < 2.1. These are the model of 
Malhotra (1993), the Nice model of Tsiganis et al. (2005) and the resonant model of Morbidelli et al. (2007). Here 
we use the configuration of Malhotra (1993), since it has the most relaxed initial conditions for the giant planets: 
the initial period ratio is PsIPj - 2.05 and the giant planets have their current inclinations and eccentricities. The 
terrestrial planets remained on their current orbits. 

Figure |2] shows the result of numerical simulations using Malhotra's planets. The layout is exactly the same as 
fig.[Tl and there are some significant differences. First the chaotic region is much narrower In the current solar system 
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the chaotic region ranges from 10° to 60°, though from to 60° on longer time scale (Laskar & Robutel, 1993; Laskar 
et al., 2004). In fig.|2] the chaotic region is confined to sq from 35° to 60°. The chaos from the secular resonance 
(gi - gs) - (ii - S2) has disappeared because ^5 > g2- The only source of chaos is the overlap of the two resonances 
ifr - \s2\ and ij/ = \si\. This can be inferred from the bottom panel of fig. |2] the lines around si and S2 are much 
narrower than in fig. [1] implying a cleaner spectrum and thus a periodic signal instead of a chaotic one. The obliquity 
can no longer resonate with all the intermediate frequencies and the chaos is confined to the main resonances. By 
increasing ^5 and breaking the resonance {g\ - g^) - {si - S2) the regions of initial obliquity with values below 35° and 
above 60° are stable, though the amplitude of oscillation is still large. This could be reduced if the terrestrial planets 
were dynamically colder, which is discussed in the next subsection. 

4.3. Compact giant planets and colder terrestrial planets 

Here we examine the stability of the obliquity of Mars in a solar system where the orbits of the terrestrial planets 
are dynamically colder 

The orbits of the terrestrial planets prior to the LHB cannot be obtained in similar fashion to those of the giant 
planets: the former were formed by solid body accretion over 100 Myr time scales well after the gas disc dissipated 
(e.g. Raymond et al. 2009). Their original orbital configuration was determined by the stochastic process of planetary 
embryo collision and accretion. Studies of the dynamical instability responsible for the LHB, likely caused by giant 
planet scattering which caused a large semi-major axis jump for Jupiter (Brasser et al., 2009; Morbidelli et al., 2010), 
cannot determine the exact effect this would have on the orbits of the terrestrial planets individually, so we study them 
as a system. One useful metric is the Angular Momentum Deficit (AMD; Laskar, 1997), which measures the deviation 
of the terrestrial system from circular and coplanar orbits, and is given by (Laskar, 1997) 



where rn„ is the mass of each planet Mercury to Mars in solar masses, the semi-major axes a„ are measured in AU 
and the inclinations ;„ are with respect to the invariable plane. Its current value is AMDq - 1.4043 x lO^^* and the 
partitioning among the planets is 35.8%, 19.9%, 20.1% and 24.2% for Mercury to Mars. This is close to the 25% one 
would expect for each planet. 

The AMD is used as the independent variable to construct a set of terrestrial planet orbits with a range of dynami- 
cal excitation, ranging from 10% to 100% of the current value (where the current AMD is equal to unity from now on) 
in ten equal steps. In order to test some low AMD values, we chose to place the giant planets on the resonant, circular 
configuration of Morbidelli et al. (2007), to remove any forcing from the giants on the terrestrials. It is assumed that 
the dynamical instability that caused the LHB would only excite the terrestrial planets, thus increasing their AMD to 
its current value. We have tested this assumption by performing simulations in which we subjected the terrestrial plan- 
ets on random orbits with pre-assigned AMD values to a so-called 'jumping Jupiter' scenario (Brasser et al., 2009). 
We consistently find an AMD increase of ~ 1 and that the increase is systematic, suggesting that the primordial AMD 
was most likely lower than its current value. This can be argued as follows. During the instability Jupiter suddenly 
appears on an orbit that is closer to the terrestrials and more eccentric, and remains there. This adds a forcing to the 
eccentricities of the terrestrials, whose amplitude depends on the eccentricity of Jupiter and the semi-major axis ratio 
between Jupiter and the terrestrials. Thus the eccentricities of the terrestrial planets, and thus the AMD, experience a 
systematic increase. Therefore we take the upper limit of the primordial AMD equal to 1 and focus our discussion on 
values closer to 0. 

For simplicity we kept the share of the AMD of each terrestrial planet the same as its current one. For the current 
orbits of the terrestrial planets, with the exception of Mercury, e ~ sin /. Therefore, the share of the AMD for each 
planet was partitioned equally among its eccentricity and inclination. The other orbital parameters remained equal to 
their current values. 



AMD = 




'■n 



Vfl^(l - ^f^^,cosi„) 



(7) 
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Figure 3: Similar to Fig.^ but now for a compact configuration of the giant planets and colder terrestrial planet orbits. The AMD for the top two 
panels is equal to 40% of the current value. The top-left panel depicts the maximum and minimum obliquity as a function of original obliquity. 
The top-right panel shows the averaged proper precession frequency. The bottom panel is a summary that shows the stable, chaotic and resonant 
regions as a function of initial AMD and initial obliquity for the systems studied here and the current solar system. 



The summary of these numerical experiments is presented in Fig. [3] taken over 256 Myr. For the top two panels 
the AMD is 40% of the current value. Note the small amplitude of oscillation of the obliquity, ~ 12°, and strictly 
confined range of obliquity motion for eo ^35° in the top-right panel. This is an indication of stable motion. The first 
resonant zone, where ij/ = \s2\, occurs for initial obliquities 37° to 45°. In resonance the obliquity oscillates between 
27° and 47° with a period of approximately 3 Myr The second resonant zone, where t}r = is confined between 
initial obliquities of 48° to 56°. Here the obliquity oscillates between 42° and 57° and with a period of 4 Myr Above 
an initial obliquity of 57° the motion is stable yet again, and the range of oscillation is merely ~ 9°. The top-right 
panel depicts ij/. There are a few very narrow chaotic zones associated with the crossing of the separatrix between res- 
onant and non-resonant motion. Inside the resonances )}risa constant and outside the resonates it varies asi}/ = a cos s. 

The bottom panel is a summary plot that shows the stable, chaotic and resonant regions as a function of initial 
AMD and initial obliquity. The situation for the current solar system has been added to the right. During the Noachian 
there are two large stable regions for mean obliquities below 35° or above 60°. In between these regions the motion 
depends on the AMD: for values less than 0.6 there are two stable resonant zones. These resonances start to overlap 
when the AMD is larger than 0.6 and the motion becomes chaotic. The chaos is confined to the resonant zones because 
there is no other source of chaos. We expect the above results to also apply when taking the giant planet configurations 
of Malhotra (1993) and Tsiganis et al. (2005). In all three models ^5 > g2 so that the resonance {g[ - gs) - {si - ^2) 
has disappeared and thus the chaos is confined to the resonant region i.e. for sq between 30° and 60°. For Malhotra's 
model, however, forcing from Jupiter limits the minimum AMD of the terrestrial planets to approximately 30% of the 
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current one, while for the other two models the minimum AMD can reach all the way to 0. 

Currently Mars' obliquity oscillates with an amplitude of approximately 20° and period 125 kyr (e.g. Ward, 1974). 
For the systems studied in Fig.[3]above, the obliquity oscillation amplitude in the lowermost stable region is also 20° 
when the AMD is equal to the current value. The oscillation amplitude then decreases as AMD'''^, so that it is just 6° 
when the AMD is 0.1. 

5. Summary and conclusions 

We have analysed the evolution of the Martian obliquity for different initial conditions of both the terrestrial plan- 
ets and the giant planets, with a focus on the compact configuration of the giant planets presented in Morbidelli et al. 
(2007). The zones of stable obliquity are summarised in the bottom panel of Fig. [3] 

The obliquity could have been chaotic if the AMD of the terrestrial planets was more than 60% of its current value, 
but only when the mean value of the obliquity ranged from 30° to 60°. For these values of the AMD, the secular spin- 
orbit resonances ip - \s2\ and ij/ - \s\ \ start to overlap. There, significant changes in the mean obliquity as well as 
large-range obliquity oscillations become possible. Outside the resonant zones the obliquity oscillates with an ampli- 
tude set mostly by Mars' share of and the total value of the AMD; the maximum value of the oscillation amplitude of 
the obliquity is approximately 20°. Before the LHB, the stable zone with obliquities lower than the resonant values 
is increased by approximately 25° when compared to the current solar system. Therefore the obliquity was stable for 
any value of the AMD as long as its mean value was below 30° or above 60°. For values of the AMD less than 0.6, 
the obliquity is stable in the resonances with 52 and si, but could suffer large-amplitude, long-period oscillations. 
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